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1. IMTRODUCTIOii 

Consider the Initial value problea for a system of ordinary dif- 
ferential equations 


y'(t) - f(t.y). y(0) - y^j. 


( 1 . 1 ) 


To solve (1*1) niBserlcally, a conventional q-stage pth-order Runge-Kutta 
(RK) method proceeds from t^ to t^^^ • t^ + h by evaluating 


ki - hf(t^ + ■*' ^ ^lj*^j-l^» * " 0* (l-2a) 

^1 


where n. * E p. , and combining the values k to yield 
■ j-1 ^ 




(1.2b) 


The coefficients 6^^ and are determined so that when the true solu- 
tion of (1.1) at t^, substituted Into (1.2) and (1.3), the 
value y should agree with the Taylor expansion of y(t ^,) at t , 

ttrl ttrl n 


2 p 

5T + ••• + ?r * 


In at least the first p + 1 terms. This pth-order accuracy Is usually 
denoted as 




(1.3) 


Solving (1.1) with a Runge-Rutta method as described above Is 
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rather expensive since It requires at least q function evaluations at 

any t . First, the k.'s are obtained to form y . , • Then, If based on 
n 1 nri 

some kind of error estimate, It Is decided that y_^. Is not accurate 

ttrl 

enough, the current stepslze h Is reduced and the k^'s are re-evaluated 
to obtain another value for y^^> This last step Is repeated until 7^]^ 
succeeds the error test* So, the number of function evaluations Is 
q + m(q-l), where n Is non-negative but not always 0* Itoreover, In 
order to keep m low, h Is usually chosen In a somewhat conservative 
manner, to a degree that the method is not as efficient as it should be. 

Fortunately, there Is a way to avoid the cost of re-evaluation. In 
developing an RK-llke technique for starting an automatic high order 
multistep ODE Integrator, Gear [3] proved the existence of values 
3.,, 1-0, . . . ,q‘*-l, j-l,.«.,l and y. , j-l,...q', s-1, ... ,p which, when 

Ij j8 

used to compute k^ from (1.2a), give pth-order accurate approximations 

8 ( s) 

of the first p scaled derivatives h y' ' (t )/sl, namely 

n 

It - £ •‘j-iTj, + • - 1 a-4) 

Since a Taylor's series method of order p takes the form 

2 p 

(1.4) naturally suggests the formulation of the following modified 
q'-stage pth-order Runge-Kutta method: 
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(1) Choose a stepslse h^. 
based on 


Calculate and the scaled derivatives 


5 + o(hfi). 


1 - 0, (1.6a) 

s*l|..»jp. (1* 6h) 


(2) Estimate *:he adjusted stepslze h based on an errc.* estimate and 
scale the derivatives In (1.6b) accordingly. 


— y^®\t ) - r® ^y^®\t ) 
si y ^ s! y 




"" ). 8-l,...,p (1.6c) 

j_l J i 38 


where r - Vh^, Vn ^ < ^max ''mln» 'max' 

to obtain y . vhlrh now can be expressed as 

ttrl 


iri-l 




E r" i: k 
8-1 j -1 


j-l^j8‘ 


(1.6d) 


y^j^ obtained In this way also satisfies (1.3), therefore Is accurate to 
a pth-order. Observe that h In (1.6c) must be determined so as to yield 
an appropriately small estimate of the 0(h^^) term. Once a reliable 
error estimator has been devised, this can be done with no substantial 
efforts since there Is no need to re-evaluate kj^. Certainly, some con- 
ditions must be lnq>osed on the ratio h/h^, but If the stepslzes are con- 
trolled properly, the modified Runge-Kutta method can be superior to a 


conventional one 
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In this thesis » questions raised by the proposed aethod are stu- 
died* The questions considered are: 

1. How to estiaate the local truncation error and accordingly adjust 
the stepslze h? 

2* What is the region of absolute stability which is now a function of 
the ratio r - What is the value of r that yields the largest 

region? 

3* How does the aodlfied Runge-Kutta aethod compare to a conventional 
one in numerical testing? 

For practical purposes^ we limit ourselves to methods of fourth order. 
Nevertheless > o could hope to get a general understanding of this 
class of methods without the necessity of going into the overwhelming 
mathematics of methods of order greater than four. 
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2. ERROR ESTIMATION 


In this section, we discuss In details as how to construct an error 
estimator required by the algorithm (1.6a-d). What we really need to 
estimate here is the local truncation error of y^^^ defined by 

p q' 

and accurate to a pth-order, that Is , 

Vl ■ 


For this purpose, we employ a technique similar to the one used In the 
Fehlberg fonmilas (see Fehlberg [2], Bettis [1].) Oiir objective Is to 
find y . which is also some combination of k. , but is of one order 
higher in accuracy. 


'nfl 


nfl 


■ "a * 


and to use y . , - y . , as an estimate of the local truncation error of 
•'nfl ■'n+1 

^n+l* 


It turns out that such exist, but the number of stages 
q" required to obtain both y^j^ and Is slightly higher than the 
number of stages needed to obtain only This is not unexpected. 
For example, the classical fourth-order RR method requires only four 
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stages while the Fehlberg formulas of order four and five need six 
stages, but the latter Include an useful error estimator for the 
fourth-order value* In Gear [4], It Is shown that fourth-order formulas 
given in (1*4) are possible with only six stages. In the following, we 
will show that with one more function evaluation, one can generate both 
fourth-order values for the scaled derivatives and a fifth-order value 
for which an error estimator can be formulated as described 
above * 


We begin by deriving fourth-order formulas using six function 
evaluations (p ■ 4, q" - 6 in (1.6a) and (1.6b), for convenience we con- 
slder f as a function of y only.) If we express the values h^y' ', 
s*l,...,4 in terms of their elementary differentials (see Gear [3], (5],) 


V' 


u2 (2) 

4 ”' 

( 4 ) 


¥• 

hjif/ + f^i. 


and expand k^, i»0, ...,5 in terms of elementary differentials of order 
up to four, for example. 


“ V’ 

kj - hj^f + a^hj^f^f + jY ^ ''' 

and so on . . • , then to obtain fourth-order approximations of the scaled 
6 f 8 ) 

derivatives h^^y' /s!, s*l,...,4, the following Identity must be 
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satisfied by and Vjgi 


A r 


1 

0 

0 

0 

0 

0 

0 

0 



0 

0 

0 

0 

6 

41 

_ 3 _ 

41 

2 

4 ! 

1 

41 


where F Is the 6x4 matrix [Yjgl A Is an 8 x 6 mat lx 

2 2 3 2 

correspond to the terms f, f2f » fj^f* f^f . ^2^1^ ’ 

respectively, and whose columns correspond to 1*0, ....S. 

to be 



( 2 . 1 ) 


whose rows 
and fjf 
A Is found 
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with Q^, leflned as 


1 

Pj ■ E 1 ■ 2«3y4|5) 

J-2 

^ 2 

Ql - E ^ - 2, 3,4, 5, (2.2) 

1 

■ 2 1 ■ 3,4,5. 


It Is clear from (2.1) that the first colvjmn of r is equal to 
T ^ 

[1 ,0,0, 0,0,0] and for s • 2,3.4, Yi " ” 2 y. • So, after eliminating 

IS ^2 

the first row and coliimn in A and the right-hand side of (2.1), there 
are only 21 conditions to be satisfied by 30 unknowns i»l,...,5, 
j«l,...,i t'nd y. , J“2,...,6, s*2,3,4. A solution of (2.1) can be 

j® 

found. In fact. Gear [4] showed the existence of a nine-parameter fam- 


ily of solutions. 


Now, in addition to the fourth-order approximations of the scaled 
derivatives, is it possible to estimate the local truncation error of 


'n+1 


4 6 

- y + E Z k 


s-lJ-1 




(2.3) 


still using the same six stages k^7 In other words, if we denote the 

principal truncation error term of (2.3) by there 

T 

exist a non-zero vector ^ such that 

6 c c 




(2.4) 
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The answer is negative and a pr^nf by contradiction Is given below. 

Assume that we can find A and T satisfying (2.1), and asaune 

further that (2.4) holds for some non-zero vector y. Since the columns 

in the right-hand side of (2.1) are linearly Independent, it follows 

that the columns of r, denoted as y' ' for sal,...,4, are linearly 

Independent. Also, since (2.4) implies that Ay « 0, y must be linearly 

Independent of y' Otherwise, there exist scalars I , not all of them 

8 

zero, such that 

E - Ay - 0, 

8-1 * 

which violates the linear Independency of toe columns In the right-hand 
side of (2.1). 

Let B be the 4x6 matrix consisting of the last four rows of A. 
Then the four Independent vectors y^^\ ^ belong to the 

null space of B. Therefore, B has rank at most 2. If - 0, the 

method reduces to a five-stage method which Is non-existent (see Gear 
[4],) so we can assume that a~ * 0. Thus the first row of B is indepen- 
dent of the remaining rows because of the zeros in Aj^2» i“6»7,8. Hence 
these rows assume rank at most 1. Since Ag^ « 0, this Is possible only 
If ®2^2 " ^2 " (2.2), we have Q2 - ‘2“l’ ® since 

* 0. Again from (2.2), this Implies that - 0, which then 

leads to * Q3 " 0 as a consequence of the linear dependency of the 
last three rows of B. For to be 0, either - 0 01c P^ • 0, but 
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not both. Otherwise > s five-stage nethod results. So assuae that 

0 and consider the following subsyswea of (2.1) 







“l 



“1 

“2 

“4 

“5 


21 

0 

0 

2 

2 

2 

2 



2 


“l 

“2 

“4 

“5 

F ■ 

0 

3! 

0 

3 

3 

3 

3 


0 


6 

“l 

“2 

*4 



0 

4! 

0 

0 

*4^ 

“5**5 


0 

0 

3 

4! 

U. 



.. 






where r Is F after crossing out the first column and the first and 

fourth rows. It Is fairly clear that the matrix In the left-hand side 

of (2.5) 1 aon-slngular. Since Af - 0, this non-singular Ity Implies 

thjL ■ 0, 1»2,3,5,6. But then, the first and fourth rows of A imply 

that ^ is identically zero, which is a contradiction. Hence we must 

have Pg ■ 0. As a result, - ^43**2 ^ ***** consequently, 

a,P, ■ Q, ■ 0. Now, P, cannot be 0. Otherwise, Is 0 and hence (2.1) 
4 4 4 4 5 

cannot be satisfied. So, 0. Again, consideration of a subsystem 

of (2.1) similar to (2.5) leads to the conclusion that ^ Is zero. Our 
proof Is thus complete. 

Since using only six function evaluations Is not enough to obtain 
the desired error estimator, we have to do some extra function evalua- 
tions, but how many? It turns out that we need only one more. When the 
coefficients given by the Fehlberg fonulas of order four and five 

are substituted In the matrix A, we discover that It Is not consistent 
with the right-hand side of (2.1), thus cannot be solved for F. How- 


II 


ever, by chaaglng the parameter In the sixth column of A, the incon- 
sistency can be reaoved. Consequently, if we still maintain a six-stage 
RK process, but evaluate the sixth stage twice: the first time to 
obtain and use it with k^, i*0, ••*,4 to solve for F; the second time 
to get k* which, together with k^, i«0, ...,4, gives rise to a vector 
y » fyj^i • • • »ygl^ so that the following 

^iri-1 * ^n *^5^6, 
holds, then our goal is achieved. 


Using this approach, a set of coefficients 3m snd y. has been 

i J J8 

calculated and is given in Table 2.1 and Table 2.2. .Vid the error esti- 
mator to be used is defined by 


5 

I k 
j-1 


i-i’j 


* Ve 


4 6 

E i: k 

8 « 1>1 


J-1 


y 


js* 


What remains to be discussed now is how to develop a step adjust- 
ment scheme which, based on the information provided by the stepslze h^^ 
and the error estimator, calculates the stepslze h that is actually 
taken. Also needed is a strategy to select h^ for the next step. At 
this time, our knowledge on the subject relies more on experimentation 
and heuristics than on conclusive mathematical analysis. So, we post- 
pone any discussion until section 4 where numerical implementation and 
testing are described. 
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Table 2«1 Coefficients 



• 

• 

1 

4 





o 


3 

9 






32 

32 





*• 

1932 

2197 

7200 

2197 

7296 

2197 



"4J 

I 

439 

216 

1 

“ 8 

3680 

513 

845 

4104 


"5J 


289 

1 4 

2072 

169 

3 


216 

3 

513 

1026 

8 

* 

^5j 


8 


3544 

1859 

11 


’27 

2 

2565 

4104 

■ M 


Table 2.2 Coefficients and y 







Y 

Y 

Y 

Y 

? 

1 

201 

1393 

137 

16 


80 

540 

144 

135 

0 

0 

0 

0 

0 

0 

7232 

114688 

3776 

6656 


1425 

12825 

855 

12825 

0 - 

15379 

81289 

10985 

28561 


4560 

10260 

2736 

56430 

A 

261 

159 

71 

9 

U 

100 

25 

20 

" 50 

0 

9 

“ 5 

24 

5 

~ 3 

2 

55 
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3. REGIONS OF ABSOLUTE STABILITY 


When a modified six-stage fourth-order Runge-Kutta method 


nfl 




4 

E 

s-1 




(3.1) 


is applied to the test equation j' ■ \y, the value at step t^^^ Is 

related to y^ of the previous step by the sli4>le e]q>ression 


y 


iH-1 


P(|A,r)y^ 


where P(p«r), considered as a function of p « hX, Is a sixth-degree 
polynomial whose coefficients are functions of r ■ llierefore, by 
keeping r fixed, we can define the region of absolute stability of (3.1) 
In the same way aa defining the region of absolute stability for a con- 
ventional RK method. It Is for a fixed r, the set of all complex values 
p for which I P(p,r) I < 1. In this section, we are Interested In plot- 
ting these regions of absolute stability for different values of r. 


To obtain P(p,r), we first make the substitution f(t,y) ■ Xy In 
For example. 


“o ■ 

^1 - + Vl > ^ 


In general, we have 
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1+1 


where 6|^^» 1-C,...,5, ■"! 1+1 are found to be 


®0« • 

1 




«1. = 

1 

*1 



«2. = 

1 

®2 

'2 


*3m • 

1 

“3 

^3 

*3 

«4. = 

1 

“4 


*4 ^ 

*5. = 

1 

“5 


Rj Vj 

with ttj^, Pj^, defined as 

In (1.2) 

and 

(2.2), and 


\ - 

1 

j-4 ^ 

-1» 

1 4,5, 


«5 - PjsV^. 


(3.2) 


Now 


t using (3.2) to replace In (3.1), we get 


4 6 j 


'■>^1 * . fi ' 

-[1+ :: I j 6 (h^x)“]y„ 

8-1 n^lj-m ■' ■* 

- P(|i,r)y„. 


Since Yj^ satisfy (2.1) and In particular, - ll»0,0,0,0,0] , It fol- 
lows that for 8-1, 


6 6 


r E E - hX - p, 

Bi-1 J— IB 
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and for s ■ 2»3»4, 


nrlj^ 

+ (''u\ * Ts,V 5KV>’ + lf6.''s<V>*l 

* vr •‘‘ ♦ + Y«.V5)r-’M= + r6.W5r'-«,‘. 


Consequently , 


n/ X J.1 2^1 3.1 4 

P(H,r) - l+ n+ ^|i + jjpU 


4 4 

+ [ ^ (Y5.V4 + rs.V )r>-5]^5 ^ , r 

8-2 8-2 


To plot the region of stability, we set 


P(p,r) - e^®, 0 c [0,2nj 


(3.3) 


and solve for p. Ihls gives the boundary of the region of stability 
since all those p for which (3*3) holds, are such that | P(n,r) | - 1. 
The boundary divides the |i-coiq>lex plane Into regions and It Is not dif- 
ficult to determine which one Is a region of stability. One can solve 
(3.3) most conveniently by finding the zeros of the polynomial 

10 

Q(p,0) - P(p, ) - e , 0 e [0,2x], 


using the Newton method for example. We first divide the Interval 

[0,2x] Into N subintervals, and then find the roots of Q(p,0^) “ 0, 

0 - nh, h - 2x/N by taking an Initial guess 

n 

^^n,(0) " Vl* 
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and iterating as followi, 

•‘n.Otfl) ■ 

Since Q is a polynomial of degree six, there are six branches of the 
boundary to be traced, each of which starts at a root of Q(p,0) * 0. We 
know that for any r, Q(0,0) « P(0,r) - 1 - 0. So, we can always start 
at the origin of the p-complex plane and plot the boundary until it 
forms a closed curve. However, In some cases, the region of stability 
is disconnected. We then must find another starting point and proceed 
as above. 

Another useful fact about the region of stability is that it is 

symmetric about the real axis. This is true because the coefficients of 

16 

P are real. So, if p satisfies (3*3), that Is P(p,r) - e , then 
P(}i,r) ■ P(|A,r) - e thus p also satisfies (3.3). There- 

fore, we need only plot the boundary In the upper half-plane. 

In Figures 3.1 through 3.10, regions of stability of (3.1) are 
drawn for various values of r ranging from 0.1 to 10000. The coeffi- 
cients used are taken from Tables 2.1 and 2.2. Of course, a practical 
range o*: r is more restricted. Still, It is Interesting to see how one 
can enlarge or shrink a region of stability just by changing r. The 
following facts are observed: 
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1. In the range of r mentioned above, there are Intervals In which the 
region of stability either expands or shrinks. Also, there are 
values of r where a whole region Is split Into separate pieces, and 
values of r where disconnected subregions are Joined back Into one. 
To best Illustrate this, various regions of stability In each 
Interval are selected and put In different figures. 

2. For r » 0.525, 0.9741 and 43, the region of stability has a large 
Intersect with the real axis, as can be seen In Figures 3.2, 3.4 
and 3.8. For r < 0.515 and 1 < r < 18, the region of stability 
consists of three disjoint regions, a large one contained In the 
negative half-plane and two smaller ones located In the positive 
half-plane. This Is shown In Figures 3.1, 3.6 and 3.7. Note tnat 
the region for r ■ 0.1 Is extremely small, compared to that of a 
conventional method. For values of r near 2, It seems that the 
regions are approaching a limiting region, but then Sifter r • 2, 
the same patt *n of expanding and shrinking recurs. 

3. As r «B, the region of stability of (3.1) becomes more and more 

Identical with that of a conventional four-stage fourth-order RK 
method. One can easily see this by looking at the expression of 
P(^,r): as r becomes large, the coefficients of and |i^ both 

approach 0. In Figure 3.10, the regions for r ■ 50, 100, 1000 and 
10000 are plotted. For larger values. It Is almost impossible for 
the eyes to distinguish between two different regions. 




from 0.1 to 0.52 
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Figure 3.2 Regions of stability for r ranging from 0.52 to 
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Figure 3.3 Regions of stability for r ranging from 0.525 to 
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Figure 3.4 Regions of stebility for r ranging from 0.7 fo 0.9741 
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Figure 3.5 Regions of stability for r ranging from 0.9741 to 
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Figure 3.6 Regions of stability for r ranging from 1.0 to 2.0 
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Figure 3.7 Regions of stability for r ranging from 2 to 34 
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Figure 3.8 Regions of stability for r ranging from 34 



Figure 3.9 Regions of stability for r ranging from 43 to 50 
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Figure 3.10 Regions of stability for r ranging from 50 to 10000 
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4. NUMERICAL IMPLEMEHTATIOM 

In addition to the error eatlaator obtained In section 2, any use- 
ful lapleaentatlon of a aodlfled slx~stage fourth-order Runge-Rutta 
method naist also Include algorlthas to adjust the stepnlze h and to 
select h^ for the next step. Due to the Halted Information we know 
about the differential equation, we are unable to offer a complete 
analysis which guarantees the best selection of h or h^. We can discuss 
Instead some practical techniques which choose h In an atteo 4 >t to con- 
trol the local truncation error. 


We recall that the local truncation error committed by taking a 
stepslze h^ can be written as 

h^4»(t„»y(V) + 0(hj) (4.1) 

where <|*(t,y) depends on the method coefficients and the elementary dif- 
ferentials of order five evaluated at (t ,y(t )). Now, suppose that the 

n n 

local truncation error of 


^ s ® 

y,., - E r® r k. 


.-I j-i 


(4-2) 


computed with the adjusted stepslze h, can be expressed similarly as 


h\(t^,y(t^)) + O(h^) (4.3) 

and suppose that the first term In both (4.1) and (4.3) Is dominant. 
Then since we have at our disposition an estimate of (4.1), we simply 
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take h so that 

5 

r II error estlnate II < e (4*4) 

for some given error tolerance e* Or, equivalently, 

h - lL[ (0 

1 1 error estimate 1 1 

where u Is a "safety” factor less than 1. Though It Is simple, (4.3) 
does not hold. The best we can hope for Is that It Is close to the true 
local truncation error. However, of^ce h Is taken, there Is no available 
test for us to know whether Is acceptable. 

Instead of (4.3), we have the following local truncation error 

h%(t^,y(V.r) + 0(h®) (4.5) 

where K^*y>i^) depends not only on the coefficients and the differen- 
tials, but also on r « l^us. If we want to use (4.4), we must 

also choose r subject to the following constraint 

< ll(Kt^,y(t^))ll* (4.6) 

The problem posed by (4.6) Is hard to solve, If not Impossible, since 
the elementary differentials are generally not. known. However, In the 
case of. linear Initial value problems, an answer to (4.6) Is as follows. 

Let j, J"l,.*.,9 denote the fifth-order elementary differentials 

*1*3*’. *2*^. *2*?**- *2**i*>^> *l*2*l*^> *1*2** *1* 


* 4 * 
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respectively. If we Include these fifth-order tens In the Taylor 
series expansion of 1-0, ...,5 and let C be the 9 x 6 matrix whose 
rows correspond to and whose columns correspond to 

1 — 0 , • « . , 5 , 


<1^2 ® 3^3 * 4^4 


C - 0 


® 2^2 ® 3^3 ® 4^4 ® 5^5 

® ®3^3 <*4^4 ffgRj I 


with a^, Q^, and defined as In the previous sections, and 


1 


i - 2, 3, 4.5, 






1 - 3,4,5. 


then we can explicitly express 4 * and 4 as 




31 


and 


.8-5 




where Y. and t . are the coaponenta of the vectora Cf and 

J 8J 

a>2,*»«,4, reapectlvely . Since (j> and ^ are ccmblnations of E_ It la 

J 

not easy to determine r ao that (4.6) la aatlafled. But, when f(t,y) la 

4 

linear In y, the E tema are all zero, except E_ . - f f. Thua, the 

5> J 5 > 9 1 

conatralnt (4.6) Is almpllfled to 


l«(p)| < |«(1)| (4.7) 

2 3 1 

where «(p) - - x^^p - x^^p - Xjgp » a cubic polynomial In p - 

The set of coefficients discovered in section 2 unfortunately yields a 

small range of solution In the neighborhood of r > 1. However, it is 

believed that coefficients can be found to enlarge this range. 


Thus, for linear problems, we choose r to satisfy (4.4) and (4.7). 
For other problems, since (4.4) Is the only piece of Information we pos- 
sess, we have no other choice but to calculate r from (4.4) with a 
rather small "safety” factor u. As for the selection of h^ for the next 
step, we can either compute a weighted average of the previous h^ and h, 
or take h^ as a scaled value of h. It should be mentioned that these 
strategies are dictated by experimentation and appear to be better than 
other alternatives that we have tried. 


We have Implemented the modified six-stage fourth-order RK method 
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and Its related techniques Into an experlaental code and done sone 
tests. The set of probleos we solved consists of the linear problem 

y' - - Xf y(0) - 1, 

Integrated over the Interval [0,1] » and the non-llnear problem 

y' - (y - sln(t)) - (y - sln(t))^ + cos(t), y(0) - 0.5, 

Integrated over the Interval [0,10]. We also solved these two problems 
using the code RKF4S, an Implementation of the Fehlberg fourth-fifth 
order Runge-Kutta method. The performances of the two codes are com- 
pared and presented in Table 4.1 and Table 4.2 for the linear and non- 
linear problems, respectively. The following statistics are shown for 
each method: 

EPS error tolerance 

NSTEP total number of steps taken 

NFCN total nimber of function evaluations 

RELERR maximum relative error for the linear problem 

ABSERR maximum absolute error for the non-llnear problem 

To maintain a certain degree of fairness In comparing the two methods, 
we used the same Initial stepslzes and the same norm as used by RKF45. 

In both problems, the modified RK method produced slightly less 
accurate results for hlgh-order tolerances. This was probably caused by 
the failure of the code to adjust the stepslze efficiently although It 
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Table 4«1 Nuaerical results for the linear problem 


EPS 


Modified 

RK 




RKF45 




NSTEP 

NFCN 

RELERR 

NSTEP 

NFCN 

RELERR 

10-1 

2 

14 

0.60 

X 

10~^ 

2 

13 

0.42 

X 

1 

o 

CM 

1 

o 

2 

14 

0.66 

X 

CM 

I 

o 

2 

13 

0.74 

X 

io“^ 

10-^ 

3 

21 

0.54 

X 

10-“ 

2 

13 

0.29 

X 

lO"® 

o 

1 

** 

4 

28 

0.50 

X 

lO""^ 

3 

19 

0.14 

X 

1 

o 

10-’ 

5 

35 

0.47 

X 

10"^ 

4 

25 

0.28 

X 

io“^ 

io-‘ 

7 

49 

0.5/ 

X 

10-“ 

6 

37 

0.30 

X 

NO 

1 

o 

10-' 

9 

63 

0.99 

X 

1 

o 

9 

55 

0.34 

X 

lO"^ 

10-» 

11 

77 

0.31 

X 

1 

o 

13 

79 

0.37 

X 

10-* 

1— • 
O 

1 

VO 

16 

112 

0.67 

X 

10-“ 

20 

121 

0.38 

X 

Oi 

1 

o 

^4 

10-1“ 

24 

168 

0.14 

X 

10”® 

31 

187 

0.39 

X 

10-1“ 










Table 4.2 Nuaerlcal results tor the non-linear problem 


EPS 


Modified 

RK 



RKF45 

NSTEP 

NFCN 

ABSERR 

NSTEP 

NFCN 

ABSERR 

10-1 

6 

42 

0.35 



8 

64 

0.24 


CM 

1 

O 

9 

63 

0.44 

X 

10"^ 

10 

76 

0.47 

X 10"^ 

10-5 

12 

84 

0.70 

X 

10"^ 

13 

99 

0.50 

X 10"^ 

10-* 

17 

119 

0.72 

X 

10”^ 

16 

107 

0.14 

X lO"^ 

lo"^ 

22 

154 

0.15 

X 

10-3 

24 

170 

0.39 

X 10"^ 

10-« 

35 

245 

0.10 

X 

10-^ 

35 

231 

0.28 

X 10”^ 

10-1 

50 

350 

0.22 

X 

10~^ 

55 

361 

0.19 

X lO"® 

10-« 

77 

539 

0.40 

X 

10-^ 

85 

546 

0.14 

X lO"^ 

10-’ 

115 

805 

0.48 

X 

io“^ 

132 

823 

0.20 

X lO"® 

10-1“ 

176 

1232 

0.69 

X 

10“S 

208 

1284 

0.18 

X lO"^ 
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took fever steps to coaplete the Integration* It Is hoped that In 
future research, sovethlng can be done to strengthen this weakness. The 
linear problea did not cause any trouble for both methods, so the modi- 
fied RK method Is more expensive to use since the cost of function 
evaluation Is 7 per step compared to 6 per step for RKF45. However, In 
the non-llnear problem where the behevlor of the solution Is not so 
predictable, the additional function evaluation paid off as RKF45 failed 
In many steps and had to retry several times. 



36 


LIST OF REFERENCES 

[1] Bettis, D* G* > Eabedded Runge-Kutta aethods of order four and five, 
Texas Institute for Coaputatlonal Mechanics Report 78-2, University 
of Texas at Austin, Texas, 1978. 

[2] Fehlberg, E., Klasslche Runge-Kutta formeln vlerter und nledrlgerer 
ordnung mit schrlttwelten-controlle und Inre anvendung auf 
warmeleltungsprobleae. Computing 6, 1970, 61-71. 

[3] Gear, C. W., Runge-Kutta starters for aultlstep methods, ACM Tran- 
saction on Mathematical Software 1980, 263-279. 

[4] Gear, C. W., Runge-Kutta starters for multlstep methods. Report 
78-938, Department of Computer Science, University of Illinois at 
Urbana-Champalgn, Illinois, 1978. 

[5] Gear, C. W., Numerical Initial Value Problems In Ordinary Dlfferen- 
tlj^l Equations , Prentice-Hall, Englewood Cliffs, New Jersey, 1971. 


